Synchronization of chaotic networks with time-delayed couphngs: An analytic study 
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Networks of nonlinear units with time-delayed couplings can synchronize to a common chaotic 
trajectory. Although the delay time may be very large, the units can synchronize completely without 
time shift. For networks of coupled Bernoulli maps, analytic results are derived for the stability of 
the chaotic synchronization manifold. 

For a single delay time, chaos synchronization is related to the spectral gap of the coupling matrix. 
For networks with multiple delay times, analytic results are obtained from the theory of polynomials. 
Finally, the analytic results are compared with networks of iterated tent maps and Lang-Kobayashi 
equations which imitate the behaviour of networks of semiconductor lasers. 

PACS numbers: 



INTRODUCTION 

Chaos synchronization is a phenomenon which is of 
fundamental scientific interest in nonlinear dynamics and 
which is being investigated in the context of secure com- 
munication and neural activity [1}^. In particular, net- 
works of nonlinear units which relax to a common chaotic 
trajectory are the focus of recent research [6l[7j. 

For many applications, the coupling between nonlinear 
dynamical units is realized by transmitting a function 
of their internal variables to their neighbours. In many 
cases, the transmission time is larger than the internal 
time scales of the units. One example are chaotic lasers 
which are coupled by their mutual laser beams |5] |8]- 
[TT] . Thus, networks of nonlinear units which are coupled 
by their time-delayed variables - including time-delayed 
self-feedback - are a subject of recent research activities 

EHEHal]. 

The theoretical investigations of chaotic networks are 
mainly based on numerical simulations. However, there 
exists a powerful method to determine the stability of the 
synchronization manifold (SM) of an arbitrary network: 
The master stability function (MSF) [29]. This method 
connects topology with function. Using the eigenvalues 
of the connection matrix, a linear equation for a single 
unit is derived which determines the stability of chaos 
synchronization for the complete network. For the case 
of time-delayed couplings, the MSF has been studied for 
few systems, only [20, 25, 32j. 

The MSF is defined as the maximal Lyapunov expo- 
nent of linear equations with time-dependent coefficients. 
In general it is not possible to derive analytic results for 
these equations because the coefficients are given by the 
chaotic trajectory of the network. Thus, one has to rely 
on numerical simulations of the linear system of differen- 
tial equations (for chaotic flows) or difference equations 
(for coupled map lattices). 

The purpose of this paper is to derive analytic results 
for chaotic networks with time-delayed couplings. Thus, 
we concentrate on a coupled map lattice, a network of 



chaotic iterated Bernoulli maps, which allows an ana- 
lytic calculation of the MSF. The corresponding linear 
difference equations have constant coefficients. There- 
fore, stability of the solution is related to the roots of 
polynomials. Results from graph theory, control theory 
and algebra will help to derive analytic statements about 
chaos synchronization. 

The linear stability equations of Bernoulli networks 
may be considered as an approximation for linear equa- 
tions with time-dependent couplings. Thus, we compare 
our analytic results with numerical results for other it- 
erated maps and with simulations of rate equations for 
coupled semiconductor lasers. In the following section 
the MSF for chaotic networks with time-delayed cou- 
plings is introduced. We concentrate on networks where 
chaos is generated by the couplings and/or self-feedback, 
similar to semiconductor lasers. Afterwards, analytic re- 
sults are derived for Bernoulli networks with single and 
multiple delay times. The last section compares these 
results with networks of iterated tent maps and of laser 
rate equations. 

MASTER STABILITY FUNCTION (MSF) 
MSF with time delay 

One powerful method to analyze the stability of syn- 
chronization in networks of coupled systems with identi- 
cal units is the master stability function (MSF) proposed 
by Pecora and Caroll [29j. A network of coupled identical 
dynamical units can be analyzed by linearizing the dy- 
namical equation around the synchronization manifold. 

In order to obtain analytic results, we restrict our in- 
vestigation to coupled map lattices, and we use an iden- 
tical function /(x), x G [0,1] for the internal dynamic, 
the self- feedback and the couplings. We extend the mas- 
ter stability function of [29j for a system without time 
delay to a system with arbitrarily many different time 
delays [27| [33H35] where each term with delay time ti 
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is weighted with a positive couphng parameter ai and a 
positive self-feedback parameter r]i . The network consists 
of units with variables xl G [0, 1], where z = 1, is 
the index of the unit and t is a discrete time step. The 
system is defined by 

M 

M N 

Without loss of generality we order the coupling terms 
with ascending delay times, so that tm is the maximum 
time delay. Each coupling delay time ri has its own cou- 
pling matrix, the normalized weighted adjacency matrix 
Gi with Y^jGi^ij = 1 and dj > [24,31]. Therefore, 
the coupling is invasive and non-diffusive, it changes the 
trajectory of the coupled system in comparison to a non- 
coupled system. We assume that the coupling matrices 
Gi commute, otherwise the MSF method cannot be ap- 
plied. The self-feedback of the system is not included in 
Gi, hence Gi^u = 0. 

Complete zero-lag synchronization x] = ... = x^ = 
St is a solution of these equations. The synchronized 
trajectory is given by 

M 

St = Vo f{st-i) + E + /(^^-^) (2) 
1=1 

The stability of the SM is determined by linearizing 
Eq. ^ in the vicinity of the SM, Eq. With 5x1 = 
xl — St we obtain 

M 

^4 = Vo f\st-i) 8x\_^ ^^Vi f\st-ri) Sxi_^^^ 

1=1 

M N 

For a network with TV nodes we obtain N coupled lin- 
ear equations with time-dependent coefficients. Since the 
coupling matrices Gi commute, we can expand the A^- 
dimensional perturbation Sxt to the common eigenvec- 
tors Wk with eigenvalues 7^^/^, k = 1, A/", / = 1, M of 
the matrices Gi. For each mode k of the perturbation, 
with 5xt = £,k,t '^fe, we obtain 

M 

&,t = E + f{St-ri)ik,t-Ti (4) 

where we have defined ctq = and tq = 1. £,k,t is the am- 
plitude of the perturbation corresponding to the eigen- 
value jk,i of the coupling matrix Gi. 

To gain analytic results we focus on a chaotic map with 
constant slope, namely the Bernoulli map which is given 



by 

/(x) = (a x) mod 1 (5) 

and is chaotic for a > 1. Since f\st) = a is constant, 
Eq. Q becomes 

M 

ik,t = ^{m^ lk,l)(^ik,t-ri (6) 
1=0 

With the ansatz £^k,t = Zk^k,o^ the whole stability prob- 
lem becomes a problem of solving the polynomial of de- 
gree TM 

M M 

with ^k,i = ivi +^/7/c,0<^- 

For each eigenvalue 7^, Eq. ^ yields tm roots which 
we label with Zk^r with r = 1, ...,tm- Our goal is to find 
coupling parameters Pk,i such that: 

1. For 7o,^ = 1, / = 1, ...,M, there exists at least one 
zo^rm, ko,r^| > 1- This guarantees a chaotic 
dynamic of the SM. 

2. For each jk {k > 1) all roots Zk^r lie inside the unit 
circle \zk,r\ < 1- This guarantees a stable SM. In 
this case, the MSF is defined as 

A = max Inlzk r\ (8) 

fe>0,r ' 

Eq. ([t]) allows an analytic investigation of chaos synchro- 
nization. Later we will discuss to which extend Eq. ^ 
is a good approximation for other iterated maps with 
time-dependent slopes and even for corresponding dif- 
ferential equations with time-dependent Jacobi matrices. 

Concept of local Lyapunov exponents 

The stability of the SM of the complete network is 
determined by the linear equation Q. This equation 
yields the Lyapunov exponents of the network parallel 
(70 = 1) and perpendicular (7/c, k > 0) to the SM. As 
we will see later, it turns out that it is useful to consider 
the contributions of individual terms in Eq. Q to the 
stability of the SM separately. Hence, we define "local 
Lyapunov exponents" to discuss these contributions. For 
example, the local Lyapunov exponent is defined as the 
maximal one of the equation 

ik,t = Vof\st-i)ik,t-i (9) 

where St is the trajectory of the complete network, in- 
cluding time-delayed terms. 
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For the case of large delay times ri one finds the follow- 
ing result: When the instantaneous Lyapunov exponent 
is positive, the network cannot synchronize [25j. Note 
that local Lyapunov exponents are not the Lyapunov ex- 
ponents of the isolated units, since the network changes 
the trajectory. Only for the Bernoulli network the linear 
equations do not depend on the trajectory, in this case 
the local Lyapunov exponents are identical to the ones of 
the corresponding isolated units. 



Graph spectrum 

The stability of the synchronization manifold (SM) is 
determined by the eigenvalues of the coupling matrices 
G|, according to the MSF Eq. (|8|. We consider matrices 
with unit row sum, Gi^ij = 1. Since we consider self- 
feedback separately, we have Gi^a = 0. 

We restrict our discussion to non-negative matrices 
p6j, Gi^ij > 0, and to completely connected graphs. 
One eigenvalue is unity, 70, ^ = 1, and according to the 
Perron- Frobenius theorem 70, ^ is not degenerate and has 
the largest modulus of all eigenvalues of Gi. Therefore, 
we order the eigenvalues such that 1 = 70, ^ > |7i,^| > 

|72,z| > ... > \lN-l,l\' 

The eigenvector for the largest eigenvalue 70 = 1 is 
wq = (1, 1). It corresponds to a perturbation parallel 
to the SM. Since we discuss only chaotic networks, the 
Lyapunov exponent Amax of the mode 70 is positive; the 
dynamics in the SM, Eq. ^ is chaotic. As we will show 
later, for large single delay time r the spectral gap A = 
1 — 171 1 determines the stability of the SM. Complete zero- 
lag synchronization is possible in the limit of weak chaos, 
Amax ^ 0, if and only if the spectral gap is nonzero, 
l7i|<l. 

The theory of nonnegative matrices [36] relates the 
eigenvalue gap to the loop structure of the corresponding 
graph: A is nonzero if and only if the greatest common 
divisor of the length of loops of the graph is unity. Hence, 
for a single large time delay, it is easy to see whether a 
network can synchronize in the limit of weak chaos. In 
fact, this result has been extended to networks with mul- 
tiple delay times ri by a self-consistent physical argument 
based on mixing information of the chaotic trajectories 
of the nodes [23j. 

For some single graphs the eigenvalues are known ana- 
lytically. Fig.[l]shows the eigenvalue gap for some graphs. 
Note that a directed triangle and a directed square have 
zero gaps, but if one connects them, the greatest com- 
mon divisor of the loops 3 and 4 is unity and the gap is 
nonzero |36|. 
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Figure 1: Schemata of bidirectionally coupled pair, spectral 
gap A = 0, bidirectionally coupled triangle, A = |, unidirec- 
tional triangle, A = 0, square, A = 0.1215 and combination 
of triangle and square, A = 0.1215. These are the values for 
equal coupling weights. 



ANALYTIC RESULTS 

In this section, we present analytic results for Bernoulli 
networks. We focus on complete zero-lag synchroniza- 
tion. According to the previous section, the stability of 
the synchronization manifold (SM) is determined by the 
roots of the following polynomials 



M 

E 

^=0 



M 



(10) 



Note that each includes we omitted the index k 
in our notation. 

The network contains M + 1 different delay times r^, 
and tm is the largest one. For each mode k of per- 
turbation we have a set of eigenvalues jk,i of the cou- 
pling matrices Gi corresponding to the delay times r^. 



Hence, we have to calculate the roots of Eq. (10) with 
A = {vi + lk,i) for each mode of perturbation where 
Gi is the coupling strength and a is the slope of the 



Bernoulli map. Eq. (10) includes the local dynamics with 
To = 1 and (To = and self-feedbacks with delay time ti 
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and strength r]i . The mode parahel to the SM has eigen- 



values 7o,^ = 1 for ah terms of Eq. (10). 



Some results can be derived immediately. The polyno- 
mial Eq. ( p^ has tm roots Zr and the theorem of Vieta 
gives Hr \^r\ = \Pm\' Hence, for |/3m| > 1 at least one 
root is outside of the unit circle. Consequently, all modes 
k with \{r]M + ctm 7k, m) > 1 are unstable. If the cou- 
pling (Jm and the self-feedback r]M of the largest delay 
time Tm are such that there exists one eigenvalue 7/c,m 
with k > and \r]M + ctm 7fc,M| > l/a, the network can- 
not synchronize, the SM is unstable. 



Eq. (10) denotes the roots of a polynomial P{z) which 
can be written as 



M 



Piz) 



(11) 



/=0 



If the mode k is unstable, P{z) has at least one root with 

\z\ > 1. Let us assume XI ^^o (^^ + 7fc,0 <^ ^ 1? which 
gives 



M 



(12) 



1=0 



Furthermore, on the real axis we have 

lim P{z) oo (13) 

Since P{z) is continuous we can conclude that P{z) has 
a root on the real axis with 2^0 > 1. Hence, the mode k 
is unstable if 



M 



(14) 



^=0 



The Gershgorin's circle theorem [37J states that all roots 
of P{z) lie inside a circle with radius R given by the 
inequality 

R < max<^ 1, ^\{m^ lk,i) ol\ > (15) 



Hence, the mode k is stable if 

M 



^\{:ni ^crnk,i)(^\ < i 



(16) 



^=0 



Eqs. (14) and (16) give the parameter regime for chaos. 



Note that this is valid for all perturbation modes k. 

For a chaotic trajectory the perturbation mode paral- 
lel to the synchronization manifold, 70, ^ = 1 has to be 
unstable. This is the case for 



M 



^{Vi ^cri)a> 1 



(17) 
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Figu re 2 : Spectrum of Lyapunov exponents for the system of 
Eq. ([20| and T = 40. a) /3o = 0.5, b) /3i = 0.5. There is at 
least one positive Lyapunov exponent for /3o + /3i > 1 



whereas for 



M 

E 

z=o 



{r]i ^ai)a <1 



(18) 



the perturbation mode for /c = is stable. Hence, the 
border to a chaotic trajectory is given by 



M 



M 



(19) 



^=0 



^=0 



If the sum of couplings is larger than 1 the system is 
chaotic. Note that this transition to chaos does not de- 
pend on the delay times ri in contrast to the region of 
synchronization which is very sensitive to the values of 
the delay times r/. In this work we consider only param- 
eters where the network is chaotic. Another conclusion 
from this result is the fact that for some networks it is al- 
ways possible to find coupling parameters and Bernoulli 
slopes for which stable chaos synchronization is possible. 
This networks have to have eigenvalues k > which ful- 
fill Eq. ([l6| for parameters which ensure 17 This is the 
case for IjkA < 1 where k > 0. 



Single delay time 

For a network with a single delay time r, including 
coupling and self-feedback we have to find the roots of 



in the chaotic region 



+ (^i + ^i7/c)ce (20) 



im^m + cri)Q^ > 1 



(21) 



^=0 



For each mode /c, Eq. (20) has r roots zj^^r which define a 
spectrum of r Lyapunov exponents Xk,r = ln|2;/e,r|- Fig. [2] 
shows this spectrum as a function of the parameters /3o 
and Pi . For Pq-^Pi > the maximal Lyapunov exponent 
is positive. Using the Schur-Cohn theorem [38j the region 
of stability can be calculated numerically. Fig. |3] shows 
the result for r = 2,4 and r ^ 00. With increasing 
delay time this region shrinks to the symmetric triangle 
< l-/3o. 

In fact, one can derive the region of stability analyti- 
cally for r ^ 00. In this limit, we have to consider two 
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Figure 3: Synchronization regime for = /3o ^""^ + /3i and 
T = 2 (light blue), r = 4 (blue), r ^ oo (dark blue). Note 
that T ^ oo regime is a subset of r = 2 and r = 4 regime and 
that T = 4 regime is a subset of r = 2 regime. 



cases: a) the Lyapunov exponent is of order 1, b) it is 
of order 1/r ^3 ^ . We omit the index k of the root 
Zk^r of the perturbation mode k since is included in 
From Eq. (I20|) we obtain 



/3i 



/3o 



(22) 



For case a) \zr\^ diverges for r ^ oo in the region of 
instability, hence we find = /3o for r ^ oo. If /3o > 
1, the SM is unstable for any perturbation mode with 
eigenvalue ^k- But /3o > 1 means that the uncoupled 
units are chaotic. Thus, we reproduce the result already 
found in [25j: If the local Lyapunov exponent defined 
in the previous section is positive, a network cannot be 
synchronized by time-delayed couplings if the delay time 
is much larger than the local time scales. 

•e^^- with = K/r. 



For case b) we write Zr = e^"^ 



Eq. (22) gives 



For r ^ oo we obtain 



(3^-2(3o cos (pr + 1 



(23) 



(24) 



For r ^ oo the phases 0^ are uniformly distributed on 
the circle [0, 27r][41] 02], hence the maximal Lyapunov 
exponent is given for 0^ = which yields the region of 
stability 

\Pi\<l-po (25) 
Consequently, the network can synchronize if for all /c > 
liVi +cri7fe)<^| <l-r]ooi (26) 



For large delay times and without self- feedback 7^1=0, 



Eq. (26) reduces to 



|cri7ia| < 1 - 7^o<^ 



(27) 



where 71 is the eigenvalue with the second largest abso- 
lute value, 1 > 1 71 1 > I72I... Neither the sign nor the 
complex phase of the eigenvalue 71 have an influence on 
the region of stability. This is different for small val- 
ues of r where the region of stability depends on the 
complex phase of 71. Note that this symmetry for large 
delay times has recently been proven for general chaotic 



networks with a continuous dynamics [30j. Eq. (27) has 
interesting consequences. The network is chaotic if 



ai a > 1 — r]o a 

It can synchronize if 

ctq a + \cri 7i a| < 1 



(28) 



(29) 



Thus, if |7i| = 1, Eq. (27) is identical to Eq. (28) and the 



chaotic network cannot synchronize for any set of param- 
eters. For bipartite networks one finds 71 = —1, hence 
bi-partite networks cannot synchronize completely, only 
sublattice or cluster synchronization is possible [24^ 31} . 
According to Fig. [T] a pair, a square without diagonals 
or any directed ring cannot synchronize completely for 
large delay times. A sublattice or cluster synchroniza- 
tion however, is possible [23l [241 [31] . On the other hand, 
a triangle with bidirectional couplings can synchronize 
completely, since 71 = — ^. 

We exemplify this result by a system consisting of two 
bidirectionally coupled Bernoulli units with one time de- 
lay r for the coupling and the same time delay r for 
self-feedback. By adding a self-feedback with delay time 
r, one can achieve complete synchronization in a bidirec- 
tionally coupled pair. The dynamical equation of unit i 
is defined as 

xl = (i-e)f{xi)+eKf{xi_,)+e{l-K)f(^xi_,) (30) 

with i,j G {1,2}. According to Eq. ( pS] ), the system is 
chaotic for any parameters f'{x) = a > 1 and e^n ^ 
[0, 1]. The pair has the eigenvalue 71 = —1 which gives 
= {I - e)a, pi = e{2K. - 1) a. Eq. ([25| yields two 
synchronization borders for r ^ 00 

- 1 



1 



a 



2ae 
1 



2a£ 



(31) 



For a triangle with corresponding equations we obtain 
with 71 = — ^ the coefficients /3o = (1 — e) a, /3i = £ a — 
^e{l — a = I (3 Ai: — 1) a which gives 



n^{e) = 1 



Hi {e) = 



(2- 



2 (a - 1) 

Sea 
-s)a-2 



Sea 



(32) 
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Figure 4: Blue dashed line: border to synchronization for 
a bidirectionally coupled Bernoulli pair with a — 1.5, self- 



feedback, 



oo and e,K notation. Thick line: , thin 



line: k~ . Red solid line: border to synchronization for a 
bidirectionally coupled triangle, same parameters as for the 
pair. 

Both results are shown in Fig. [4] Without self-feedback 
= 0), a bidirectionally coupled pair cannot synchro- 
nize whereas a triangle synchronizes for £ > 2 (a — 1) /a. 

Note that the upper boundaries do not depend on 
the delay time r, whereas with decreasing delay time the 
lower boundary moves down, the region of synchroniza- 
tion increases. 



From Eq. (24) a relation between the MSF and the 



maximal Lyapunov exponent in the limit of r ^ oo can 
be derived. We consider a network without self-feedback 
and with local stability /3o < 1. The MSF, i. e. the largest 
transversal Lyapunov exponent is given by 



A 



1^1^71 1 



l-/3o 



(33) 



whereas the maximal Lyapunov exponent, which is the 
Lyapunov exponent parallel to the SM, is 



Hence, one obtains 

A = A 

Thus, the SM is stable if 



1 |(Jia| 



max + -ln|7l| 
T 



l7i| < e" 



(34) 



(35) 



(36) 



This equation has a fundamental meaning. For any net- 
work with a stochastic coupling matrix G, it relates the 
eigenvalue gap 1 — |7i| to the synchronizability of the 
network. If the second largest eigenvalue is smaller than 
the largest one, |7i| < 1, the SM is stable for sufficiently 
weak chaos inside of the SM, i. e. in the limit Amax 0. 



We believe that Eq. (36) holds for any chaotic net- 



for corresponding delay differential equations. In fact. 



Eq. ( 36 ) has been derived for networks with periodic dy- 
namics [43j and our numerical results of laser equations 
confirm this condition on chaos synchronization (see the 
following section). 

Eq. (36) has been derived in the limit r ^ oo. For 
finite r one either has to solve Eq. (20) numerically or 



one can calculate regions of stability with the Schur-Cohn 
theorem [38j. But even in this case of finite r an ana- 
lytic result is possible. The border of stability is given 
by = 1 where \zr \ is the maximal root for the pertur- 
bation modes k transversal to the SM. With Zr = e^"^^ 



one can decompose Eq. ( 20 ) into the real and imaginary 
parts 

cos(0^ r) = /3o cos[^^ (r - 1)] + /3i 

sin((/)^ r) = /3o sin[(/)^ (r — 1)] (37) 



The solution of Eqs. ( [37| ) gives boundaries /5i(^) and 
/3o ((/)). The maximal value of /3o(0), i.e. the tip of the 
phase diagram of Fig. [3} is obtained in the limit ^ 0, 
which gives 



1 = /?0 + /?1 

r = /?o (r - 1) 



(38) 
(39) 



The first part agrees with Eq. (19), the second gives the 
result 



max /3o 



r- 1 



(40) 



In the limit r ^ oo one can synchronize only for /3o < 1, 
i. e. if the local Lyapunov exponent (see previous section) 
is negative. But for finite delay times even chaotic units 
(/3o > 1) can synchronize if 

/3o < ^ (41) 



Two delay times 

Now consider a network with two delay times t\ and 
^2 {j2 > ^\)' Each delay time belongs to a coupling 
matrix G\ and which have identical eigenvectors, i. e. 
identical modes of perturbation. The stability of the SM 
is determined by the polynomials 



■/?2 



(42) 



We consider only networks where /3i and ^2 belong either 
to a coupling or a self-feedback, therefore we rule out 
r]i ^ and cr^ ^ for the same system. Each mode 
of perturbation has an eigenvalue 'yk,i of the coupling 
matrices of Eq. ([T]), which gives three possibilities for 



work of the structure defined in the previous section, even 



Eq. (42): a system with self- feedback ri and coupling 
r2, a system with self- feedback r2 and coupling ri and 
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Figure 5: Synchronization regime for z^'^ — /3i z^"^ + /32 
and r2/ri = p/q with p, ^ relatively prime and a) | = 2, b) 
I = 3, c) ^ = |, d) I = 4. Region (I) denotes the region 
parallel to t^e synchronization manifold, the other regions lie 
transversal to the synchronization manifold. 



a system with two couplings ri and r2. For these three 
cases we obtain the equations 





= r]oa z^'^ ^ - 


^r]iaz^^-^' - 




(43) 




= rjoaz^^''^ - 


VcJia-fk^iz^^ 




(44) 


z^' 


= r]oa z^'^~^ - 






(45) 



For simplicity, we start the discussion omitting the lo- 
cal term, 770 = 0, and considering a pair of coupled units 



with 7o = 1, 7i = — 1. Hence, Eq. (|42|) is reduced to 

- f3i z^'-^' 



z 



■P2 



(46) 



Let ji be the greatest common divisor of T2 and ri. We 
can substitute w = z^ io obtain 



(47) 



with T2 = P Ti = qji and (p, q) are relatively prime. 
Now, all roots Zr for A: > lie inside the unit circle if 
and only if all roots Wr for /c > lie in the unit circle, 
hence only the ratio T2/T1 = p/q determines the stabil- 
ity of Eq. (46). Fig. [5] shows the regions of stability of 
Eq. (46 ) for different values of p and q calculated with the 



Schur-Cohn theorem [38 . As shown before, the system 
is chaotic if (r^i + ai) a + (7^2 + cr2) <^ > 1- For the mode 
70 = 1 this means /3i + /32 > 1. For example, the point 
in region (I) in Fig. [5|i) belongs to a chaotic system. 

The stability of the SM for a bipartite network, for ex- 
ample a pair of units, is determined by the mode 71 = —1, 



which changes the sign of /3i and/or ^2 of the point in 
region (I), depending whether this term belongs to a cou- 
pling or to a self-feedback. For self-feedbacks one has 
Pi = r]ia whereas the coupling gives Pi = —aia. For ex- 
ample if Ti belongs to a self-feedback, /3i is positive, since 
T] and a are positive. For r2 belonging to a bipartite cou- 
pling, G2 has an eigenvalue 71 = — 1, so /32 is flipped 
to —P2 which gives the point in region (IV) in Fig. |§i) 
which is stable. When both delay times ri and r2 belong 
to the coupling, Gi and G2 are coupling matrices with 
71 = —1 and we obtain the point in region (III) which is 
stable, as well. In both cases we find complete zero-lag 
synchronization. But when r2 belongs to a self-feedback 
and Ti to a coupling we obtain the point in region (II) 
which is unstable. In this case we can never achieve syn- 
chronization for any parameters <Ji and <J2 for which the 
system is chaotic. 

Obviously, the symmetries seen in Fig. |5] prevent syn- 
chronization for some cases depending on the values of 
p and q. These symmetries of Eq. ( [42| ) can be derived 
as follows. If we change the sign of the roots, w = —v 



the conditions of stability do not change, but Eq. (47) is 
changed to the polynomial 



i3i{-iyvP-^^(32{-ir 



(48) 



Now, we have to consider three cases: 

If q is even and p is odd, the phase diagram has the 
reflection symmetry P2 ^ — /^2, see Fig. [sj)) for p = 3 
and q = 2. Thus, if ri belongs to the self- feedback and 
r2 to the coupling, the system cannot synchronize, since 
the condition for stability of the SM is identical to the 
condition for chaos. 

If q is odd and p is even, the phase diagram has the 
reflection symmetry Pi — see Fig. [5^) and [sji). 
Thus, if Ti belongs to the coupling and r2 to the self- 
feedback, the system cannot synchronize. 

If p as well as q are odd, the phase diagram has the 
point symmetry Pi —Pi and P2 —P2' Thus, if 
both delay times belong to couplings, the system cannot 
synchronize. 

Therefore, the symmetries of the roots of the polyno- 



mials Eq. (46) rule out some ratios of the two delay times 
for which synchronization can occur. For the ratios which 
are not forbidden by symmetries, synchronization is pos- 
sible in a limited paramter region, shown by the dark 
regions of Fig. |5] for a pair of Bernoulli units. From nu- 
merical calculations, we observe that this region shrinks 
to zero when p and q increase and we will later show that 
zero-lag synchronization is not possible if ri and r2 are 
large with a small difference. 

These results are in agreement with [SF, '44] where 
a pair with multiple feedback and multiple couplings 
with different delay times was analyzed. The time de- 
lays which lead to zero lag synchronization follow [35l 



Eq. m],Zfi\liNd^ 



0, where Nd^ are 
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Figure 6: Synchronization regime for z^"^ — /So z'^'^~'^ -\- 
/3i + and rs/n 2 with f3o 0.3. It is a) T2 4, 

b) T2 = 8, c) T2 = 16, d) T2 = 100. The symmetry found for 
T2/T1 = 2 and /3o = is restored for /3o > and large values 

of Ti. 



the delay times of the Mg different self-feedbacks, Nc- 
are the delay times of the Mm different couplings, and 
/i, ruj are whole numbers with a restricted set of possible 
values which are specific for each system. 

Note, however, that Fig. [5] holds for any network with 
eigenvalues 7^. For example, if ri and r2 belong to the 
mutual couplings of a triangle without self-feedback we 
have 7o = 1 and 71^2 = —1/2. Hence, the point in 
region (I) of Fig. [5^) is mapped to Pi — /5i/2 and 
— /^2/2. This means that the triangle can sychro- 
nize for any ratio p, q since the perturbation modes are 
located in the interior square which is always stable. This 



result is due to Gershgorin's circle theorem (15) and the 



conclusion found earlier in this section. Symmetries rule 
out synchronization for bipartite networks, only. 

Up to now we have neglected the local term 7^0 az^'^~^ 



in Eq. (42). We will show that the symmetries still hold 
in the limit r2 ^ 00, ri ^ 00, r2/ri = p/q. The stability 
of the SM is determined by the roots of 



r]oaz' 



-Pi 



(49) 



The border to synchronization is necessarily given by 
|2:| = 1 where \z\ is the maximum of all roots z^. With 
z = e^ we obtain 



l = 7?oc^e-^^ + /3i e 



-i 0r2 qlv . 



-l<pT2 



(50) 



In the limit of r2 ^ 00 we use the fact that the phases (j) 
of the roots are uniformly distributed in [0, 27r] [4T( [42]. 
Hence, there exists one root which is close to (j) = nn— = 



TiTT^ with some integer n. Following this root in the limit 
r2 ^ 00 we obtain 

l = 7/oa + /3ie-^"-^ + /32e-^""^ (51) 

Now the border to chaos 



1 = a + (7/1 + cTi) a + (7/2 + 0-2) a 



(52) 



can be mapped to Eq. (51 ), depending on whether ri and 



r2 belong to self- feedback or coupling delays. For exam- 



ple, if Ti and r2 belong to a coupling, Eq. (51) becomes 

1 = ao a - ai ae-^^^^ - cr2 ae"^^^^ (53) 

If both p and q are odd, we choose e~'^^^ = — 1 and 
Eq. (52). As synchronization regions calculated with 



the Schur-Cohn theorem show, Eq. (51) is the border to 



synchronization. Hence, synchronization is not possible. 
This shows that the ratios r2/ri = p/q for which syn- 
chronization is ruled out by symmetry do not depend on 
the local term in the limit of large delay times ^ 00. 
In fact, the numerical simulations Fig. [6] show that the 
symmetries of the complete phase diagram are the same 
as the ones proven before for 7^0 = 0, although the phase 
diagram depends on the strength 7^0 of the local term. 
We have shown that the symmetries of the stability 



equation, Eq. (49), rule out some ratios of the delay 



times. In fact, these results support self-consistent ar- 
guments for general chaotic networks. These arguments 
are based on the fact, that the information of each tra- 
jectory of each unit has to mix after multiples of time 
intervals r in order to achieve zero-lag synchronization 
[23J. 

Finally, we consider the question to which extend the 
network is sensitive to detuning the delay times r2 and ri . 
It turns out, that synchronization is extremly sensitive to 
a tiny detuning for large delay times. For simplicity, con- 
sider the case 7/0 = 0, ri = r, r2 = r + A with r ^ 00 
where A remains finite, ri belongs to a self-feedback, 
and r2 to the coupling of a bipartite network, for exam- 
ple a pair of Bernoulli units. The boundary to chaos is 
determined by 



1 < r]i a -\- a2 a 



(54) 



and the stabilty of the synchronization manifold is given 
by the polynomial for 71 = —1 



r]iaz 



(72 a 



(55) 



The stability is determined by the largest root z^ 
e^/'^ which gives 



^2A 



7^1 a e 



!) A 



Cr2 a 



= (7^1 a cos A — (72 a) ^rjla^ sin^ (j) A 

= riia ^ (72 OL — 2 7^1 (72 a cos (j) A (56) 
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Figure 7: Cross-correlation for a bidirectionally coupled 
Bernoulli pair with self-feedback delay time n and coupling 
delay time T2 — 10000. A = T2 — n. The slope a = 1.1. 



In the limit r ^ oo the distribution of the phase (j) of 
the roots is uniform in [0,27r], hence we always find one 
root with cos((/) A) ^ — 1 which according to ([54| gives a 
positive transverse Lyapunov exponent A/r. Thus, even 
when the self-feedback delay differs from the coupling 
delay by a single time step, r2 = n ±1, a pair of Bernoulli 
units cannot synchronize for r2 ^ oo. In fact. Fig. [7| 
shows that the cross-correlation immediately decreases 
to zero when the delay times are detuned. 



i. e. we discuss the limit tm ^ oo with finite r/, / ^ M. 



Eq. (57) can be rewritten as 



M-l 



(58) 



^=0 



If this equation has a root \z\ > 1 in the limit of tm oc, 
this root is obviously determined by the first term of 



Eq. (58). Hence, if the mode k is unstable for the net- 



work without long delay, the last term of Eq. (58) cannot 



stabilize this perturbation. If, however, the largest root 



of Eq. (58) approaches the unit circle for tm oo, the 



long delay has an influence. We rewrite Eq. (57) as 



(59) 



In the limit of tm ^ oo, we make the ansatz z 
q^/tm gi which gives 



Multiple delay times 

Some of the results of the previous paragraphs can im- 
mediately be extended to a network with M delay times 
[35, 44J. In particular, we can rule out complete zero-lag 
synchronization for some networks in the limit of large 
delay times. 

First: If one delay time is much larger than all other 
ones, we find the following result: If the network without 
the long delay does not synchronize, this network cannot 
be synchronized by adding the long delay. 

Second: If for a pair of coupled units all delay times 
are much larger than the time scales of the isolated units, 
then synchronization is ruled out for specific ratios of the 
delay times. 

Note that we restrict our discussion to networks de- 
fined in the previous section: To each delay time ri there 
exists a coupling matrix Gi with constant row sum, and 
all matrices Gi have identical eigenvectors with eigenval- 
ues -fk,i; k = 0, ...,A^ - 1; / = 1,...,M. The stability of 
the SM is determined by the polynomials 



M 



1=0 



(57) 



Long delay time tm 



First we consider a network where the largest 
delay time tm is much larger than all other delays r^. 



\PM\ 

^=0 



(60) 



As before, the Lyapunov exponents of the /c-mode de- 
pend only on the modulus of the coupling — 
{j]M + ctm 7/c,m) However, the maximal Lyapunov ex- 
ponent is not given by = 0. Only if all parameters 
are positive, we obtain the boundary to stability as 



\fi- 



M\ 



M-l 



^=0 



(61) 



In this case all /3i belong to self- feedbacks, only /3m is a 
coupling. 

The stability of the SM is determined by negative or 
complex eigenvalues jk^i- For this case we find an inter- 
ference of different phases (jyri in Eq. (60) resulting in a 
non-obvious value for the maximal Lyapunov expo- 
nent Amax- For example, for ri = 3 and r2 = 300 we ob- 
tain the phase diagram of Fig. [Sj This means, that a pair 
of units which is chaotic (cross in Fig. [8| can be synchro- 
nized if Ti belongs to a coupling and r2 to a coupling or a 
self-feedback. If, however, ri belongs to a self-feedback, 
r2 cannot synchronize this pair. Note that the situation 
is different for networks with |7i| < 1 , e.g. triangles. 
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COMPARISON WITH OTHER SYSTEMS 




-I.5L 



Figure 8: Synchronization area of a system with two time 
delays n = 3, r2 = 300 and fio = 0.2. The red hnes mark the 
hnes 1 — /3o = + |/32|. The black crosses exemplify a pair 
of units which is chaotic and can synchronize if n belongs to 
a coupling and r2 to a coupling or a self-feedback, but not if 
n belongs to a self- feedback. 



Symmetry 

Similarly to the case of two delay times discussed 
before, the symmetry of the polynomial Eq. (57) rules 
out synchronization for a pair of units with specific 
ratios of the delay times ri. This symmetry holds for 
= for general values of r/, but for /3o 7^ only if all 
values of ri are large. For simplicity, we discuss the case 
/3o=0. 

Firstly, we consider the greatest common divisor r of 
all delay times r/, i.e. we define pi = ri/r where the 
integers pi are relatively prime. Secondly, we substitute 



w = in Eq. (57) and obtain 



M 



(62) 



^=0 



If a root Zr of Eq. (57) lies inside resp. outside the unit 



circle, the root Wt of Eq. (62) lies inside resp. outside as 
well. Hence, the stability of the perturbation k can be 



discussed using Eq. (62). 



The system is chaotic, Eq. (62) has at least one root 
\zr\ > 1 for 7o = 1, i.e. for 



M 

E 

^=0 



{r]i ^ai)aw^ 



(63) 



If we restrict our discussion to a pair of units where each 
Ti either belongs to a self-feedback /3i = a or a cou- 
pling /3i = —ai a, we can map the stability of the SM to 



Eq. ( 63 ) for specific values of pi . We find the following 



result: Synchronization is ruled out if pi is odd for a cou- 
pling and if pi is even for a self-feedback. These results 
have been previously seen in ^35^ .44j . 



The analytic results of the previous section were ob- 
tained for networks of Bernoulli units. For this case, the 
linear equations describing the stability of the synchro- 
nization manifold (SM) have constant coefficients. This 
fact allows a stability analysis with polynomials of degree 
tm, the largest delay time. 

For networks of general nonlinear units, however, the 
stability equations have time-dependent coefficients gen- 
erated by the chaotic trajectory of the SM. For example, 
for general coupled map lattices the stability of the SM 
is determined by Eq. (|4|, where the time-dependent co- 
efficients are given by the dynamics on the SM, Eq. ([2]). 
To our knowledge, these equations cannot be solved an- 
alytically. Hence, we compare our analytic calculations 
of the section "Master stability function" with numerical 
simulations of the skewed tent map 



/: [0,1] ^[0,1], fix) 



l-x 
1-a 



if X ^ a 
if X > a 



(64) 



with < a < 1. We choose a = 0.86 which results in the 
same Lyapunov exponent for an isolated unit as the one 
of the Bernoulli map 



/: [0,1] ^[0,1], fix) 



mod 1 



(65) 



Coupled map lattices are special dynamical systems. 
Thus, it will be interesting to compare the analytic re- 
sults of the section "Master stability function" with sys- 
tems of nonlinear differential equations. In particular, 
we consider the Lang-Kobayashi (LK) equations for cou- 
pled semiconductor lasers [23| |45] . In some cases we even 
compare our results with experiments on chaotic semi- 
conductor lasers. 

For networks of differential equations, the mathemati- 
cal structure corresponding to Eq. ([T]) is defined by 

M 

f,(t) = F[x,{t)] H[xi{t - n)] + 

1=1 

M N 

Now x{t) is a multidimensional vector. For example, for 
the LK equations x{t) contains the real and imaginary 
parts of the envelope of the electric field and the popula- 
tion inversion of the charge carriers. Details are given in 
the appendix. The dynamics of the SM is given by 



M 



m = F[s{t)\+Y.^m+'yi)H[s{t- 



(67) 



The stability of the SM is described by linear equations 
for each mode k with eigenvalues ^k,i of the coupling 
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matrix Gi. The equations corresponding to Eq. Q are 



M 

+ ^ (t?^ + (Ji 7fc,0 D^[^(^ - n)] ik{t - ri) (68) 

DF and DH are the Jacobian matrices of F and H 
respectively evaluated at the SM. The master stability 
function for the LK equations is defined in the appendix. 
In the following paragraphs we consider networks with a 
single delay time, M = 1, and with double delays, M = 2. 



Networks with a single delay time 

In the previous section, we obtained the analytic result 
that the eigenvalue gap of the coupling matrix determines 
the stability of the SM in the limit of large delay times 
r, which will be realized in this section by r = 100 ns 
for the LK equations and r = 100 for the tent map and 
the Bernoulli map. The relation |7i| < exp(— Amax^) 



of Eq. (36) is the condition for a stable SM; it relates 
the maximal Lyapunov exponent Amax of the SM to the 
second largest modulus of the eigenvalues of the coupling 
matrix G. 

As a consequence of Eq. (36), a pair of units with- 
out self- feedback cannot be synchronized, since 71 = —1. 
This result agrees with experiments on semiconductor 
lasers [46j. Two lasers coupled by their mutual laser 
beams cannot be synchronized with zero time lag, only 
high correlations with a time shift of r have been ob- 
served. 

With self- feedback, however, the situation is different. 



Eqs. (26) and (31) show that a pair of units can be syn- 
chronized if the local Lyapunov exponent is negative. 
Again, this result agrees with experiments on semicon- 
ductor lasers where the self-feedback is realized by exter- 
nal mirrors [47H52] . 



Thus, qualitatively, the relation Eq. (36) is in agree- 
ment with experiments on lasers. Unfortunately, up to 
now, experiments on larger networks of coupled lasers are 
not reported. Hence, the quantitative comparison with 
lasers has to rely on numerical simulations of the LK 
equations. 

Fig. [9] shows the numerically calculated maximal Lya- 
punov exponent A^ax as a function of the coupling pa- 
rameter a for a network of lasers and as a function of 
the coupling parameter e for a network of tent maps and 
Bernoulli maps. Self- feedback is suppressed (7^ = 0). The 
parameters for the LK equations are defined in the ap- 
pendix. 

For a triangle with bidirectional couplings, it is 71 = 
Eq. iSm predicts a transition to synchronization for 
n(^) ^ 0.69, the horizontal dashed line in 



_ 1 
2 • 

Amax 'T 





Figure 9: Maximal Lyapunov exponents of the synchroniza- 
tion manifold for a network of (a) lasers modelled by Lang- 
Kobayashi equations and (b) tent maps (dotdashed line) and 
Bernoulli maps (solid line) 
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Figure 10: Zoom of Fig. |9] the critical coupling strengths 
(dashed lines) according to Eq. (|36]), for the synchronization of 
a triangle of (a) lasers modelled by Lang-Kobayashi equations 
and (b) tent maps (dotdashed line) and Bernoulli maps (solid 
line) without self-feedback 



coupling £c ~ 0.9 for the tent maps, vertical dashed line 
in Fig. 10 3), and two critical couplings cjc^i ^ 0.45 ns~^ 
and (Jc^2 ^ 45ns~^ for the three lasers, vertical dashed 

shows the cross-correlations ob- 



line in Fig. 10 i). Fig. 
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tained from numerical simulations of the corresponding 
triangle. We find the measured critical couplings of the 
tent maps and the LK equations to be in good agreement 



with the predictions of Eq. (36). 



A more challenging test of the condition Eq. ( 36 ) is a 
network with directed couplings, 
value 7i is a complex number. 



In this case the eigen- 
For example, for the 
square with one diagonal of Fig. [T] one can add an ad- 
ditional coupling strength p for the diagonal and obtain 
the eigenvalue gap of Fig. 12 The largest gap is obtained 

-IHbi Vll 



Fig. [Tok). The numerical results of Fig. 10 give a critical 



for p = I with 7i 

Combining this result with Eq. (36) and Amax(cr) 
(Fig. [9|i), we obtain the phase diagram of Fig. 13 for 
the corresponding laser network. Fig. [13] also shows pa- 
rameters for which complete zero-lag synchronization is 
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Figure 11: Critical coupling strengths (dashed lines) for syn- 
chronization measured in numerical simulations of a triangle 
of (a) lasers modelled by Lang-Kobayashi equations and (b) 
tent maps (dotdashed line) and Bernoulli maps (solid line) 
without self-feedback 



In 




Figure 12: Eigenvalue |7i|(p) for the square with directed 
couplings and one diagonal with coupling strength p, dashed 
line: position of the largest eigenvalue gap 



achieved in numerical simulations of the complete laser 
network. We define complete zero-lag synchronization 
between the lasers for isochronal cross-correlations larger 
than 0.99 in between the power drop-outs of the low fre- 
quency fluctuations which happen on a time scale of the 
order of magnitude of 10 r. While it makes sense to 
speak of complete synchronization for such values of the 
cross-correlations, the maximal Lyapunov exponent is 
still slightly positive for correlations around 0.99. Hence, 
some of the shown points of complete synchronization in 



Fig. [13] lie outside of the predicted stability border which 
corresponds to a maximal Lyapunov exponent of exactly 
zero. The quantitative agreement between the relation 



Figure 13: Phase diagram for a square of lasers with directed 
couplings and one diagonal with coupling strength p, crosses: 
parameters for which complete zero-lag synchronization is 
achieved in numerical simulations (only synchronization for 
chaotic dynamics on the SM is shown), dots: prediction of 
the stability border for synchronization, dashed line: position 
of the largest eigenvalue gap 



Networks with two delay times 

A pair of units without self-feedback cannot be syn- 
chronized if its coupling has a single delay time. In the 
previous section, however, we have shown in agreement 
with [35l Hi] that zero-lag synchronization is possible if 
the coupling contains two delay times ri and r2. Only 
if the ratio ^ is a ratio of odd relatively prime integers, 
^ = -, synchronization is excluded. The parameter re- 
gion of synchronization is largest for small values of p and 
(1- 

This result is in agreement with recent experiments on 
semiconductor lasers [53]. For ^ = f complete synchro- 
nization was observed. Cross-correlations were large for 
^e{|,|,|,|}andlowfor ^G{i,|,f}. 

The analytic result was based on the symmetry of the 
phase diagrams. Figs. [5] and [6j Here we show that these 
symmetries can be observed for the MSF of laser net- 
works and tent maps as well. With two delays, Eq. (67) 
and (68) reduce to 



Eq. (36) and the phase diagram Fig. 13 of complete zero- 



lag synchronization is remarkable. 



k{t) = F[s{t)] + (Ji H[s{t - n)] + (72 H[s{t - rs)] (69) 

+ (71 7/e,i BH[s{t - n)] ^k{t - Ti) + 

+ (T2 7fc,2 BH[s{t - rs)] ^k{t - T2) (70) 

Because of the invasive nature of the coupling, we obtain 
a different trajectory s{t) on the SM for each pair of 
coupling strengths ((71,(72) which in each case gives a 
different linear stability equation. In order to compare 
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Figure 14: Region of stability (dark gray regime) for a pair 
of (a) lasers modelled by the Lang-Kobayashi equations and 
(b) tent maps in comparison with a pair of Bernoulli maps 
(light gray regime) with two coupling delay times n and T2 
and = f , dashed line: reflection of the top-right quadrant's 
stability border into the other quadrants 
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Figure 15: Sensitivity of cross-correlations to detuning of the 
delay times T2 and n = T2 + A for a pair of Lang-Kobayashi 
equations. 



SUMMARY 



with the previous section, we fix cji and (72, vary ^^,1 
and 7/c,2 and calculate the maximal Lyapunov exponent 
of the linear equation Eq. (70). In this interpretation of 
Pi = en 7fe,i and j32 = cr2 7k,2^ the results are universal 
in view of the fact that they make a statement about the 
stability of all modes of every possible network with real 
jk,i and 7/e,2 for the chosen coupling strengths ((11,(72). 
The point 70,1 = 70,2 = 1 is shared among all networks; 
it determines the stability of the dynamics on the SM. 

Fig. [T4I sh ows the results of the numerical simulations 
of Eqs. (|69|) and ([70| with n = 20ns,r2 = 40ns,(7i = 
(72 = 45ns~^ for the LK equations and ri = 200, r2 = 
400, e = 0.9, = 0.5 for the tent map and the Bernoulli 
map. As shown in the previous section, the region of sta- 
bility should have reflection symmetry at the horizontal 



axis. The stability region for laser networks (Fig. 14 i) 
is in agreement with this symmetry. The correspond- 
ing stability regions of the networks of tent maps and 
Bernoulli maps (Fig. [Mjp) coupled with two delay times, 
show this reflection symmetry as well. 

In the previous section, we showed that synchroniza- 
tion is sensitive to detuning of the ratio of the delay times, 
see Fig. [7| Fig. [15] shows that synchronization of two 
lasers is destroyed if ri and r2 differ by about 10 ps, which 
corresponds to the coherence length of the chaotic lasers. 
The coupling has a fixed delay time of r2 = 100 ns and a 
strength of (72 = 20ns~^. The self- feedbacks have a delay 
time of Ti = r2 + A and a strength of (72 = 30 ns~^. Thus, 
in agreement with the analytic results for Bernoulli net- 
works, lasers are sensitive to detuning of the delay times 
as well dHl [49l [54] . Of course, the detailed structure of 
the cross-correlations of Fig. [15] depends on the details of 
the laser dynamics which cannot be predicted by iterated 
maps. 



Chaos synchronization of networks of identical non- 
linear units with time-delayed couplings is investigated. 
Although the units are coupled with long delay times, 
they synchronize to a common chaotic trajectory with- 
out time shift. For rather general networks with multiple 
delay times the method of the master stability function 
allows to relate the stability of the synchronization mani- 
fold (SM) to the eigenvalue gaps of the coupling matrices. 

For networks of iterated Bernoulli maps the stabil- 
ity of the SM is calculated analytically in the limit of 
large delay times. The theory of polynomials allows to 
calculate phase diagrams of chaos synchronization and 
to derive their symmetries. Finally, these analytic re- 
sults are compared with numerical simulations of iterated 
tent maps and rate equations for semiconductor lasers 
(Lang-Kobayashi (LK) equations). Some results can even 
be compared with recent experiments on semiconductor 
lasers. 



For a single delay time, Eq. (36) is the most impor- 



tant result. It relates the eigenvalue gap of the coupling 
matrix of the whole network to the Lyapunov exponent 
of the trajectory of a single unit with feedback. It is ex- 
act for Bernoulli networks with a long delay time. But 
it compares well with our numerical results for networks 
of tent maps and LK units, too. Even the phase di- 
agram of a directed network with complex eigenvalues, 
which has been calculated for lasers, is in good agree- 
ment with this fundamental equation. For a pair of units 
without self-feedback and for any bipartite network the 
eigenvalue gap is zero, hence these networks cannot syn- 
chronize completely in the limit of a large delay time, 
only cluster synchronization is possible. 

For networks with several delay times, we could not 
find a simple relation for the stability of the SM. How- 
ever, the theory of polynomials showed some symmetries 
of the phase diagrams for Bernoulli networks. Our nu- 
merical results for tent maps and LK equations showed 
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these symmetries, as well. For a pair of units coupled 
by multiple delay times these symmetries have interest- 
ing consequences. Only for special ratios of the delay 
times synchronization is possible in the limit of long de- 
lay times, in agreement with self-consistent arguments 
related to mixing of information [23j. Again, this ana- 
lytic result is in agreement with numerical simulations of 
the LK equations and even with experiments on semicon- 
ductor lasers [53j. 

These results show that networks of iterated Bernoulli 
units have universal properties. On the one hand, we 
have analytic tools to calculate the stability of chaos 
synchronization since the linearized equations do not 
contain the chaotic trajectory. On the other hand, we 
have either numerical simulations of the linearized differ- 
ence/differential equations containing the chaotic trajec- 
tory as input, or we have direct simulations of the com- 
plete network. We have found that there is good agree- 
ment between these different systems, sometimes even on 
a quantitative level. 



The Lang-Kobayashi equations and their master 
stability equations 



The Lang-Kobayashi equations in their complex form 
are 



+ ai r(t-Ti)e-'"«"i + 

N 

+ (72 Y.Gij£^{t-T2)e- 

i=i 

n*(t) = (p-l)Jth-7n'(i)- 



[T + G^n\t)\ \£\t)[ 



(71) 



(72) 



where £^(t) is the envelope of the complex electric field 
and n^(t) is the renormalized population inversion of the 
charge carriers of laser i. The used constants are listed 
in Tab. [H 

For the numerical simulation of the equations, we make 
the ansatz 



(73) 



and thus obtain the real- valued differential equation sys- 
tem 



Table I: Used constants in the simulation of the Lang- 
Kobayashi equations. Values are taken from |55] . 



Linewidth enhancement factor 
Differential optical gain 
Laser frequency 
Pump current relative to Jth 



a 5 

Gn 2.142 X 10^ s-^ 



27r c 



635 nm 

p 1.02 

Threshold pump current of solitary laser Jth 7 A^soi 
Carrier decay rate 7 0.909 x 10^ s~^ 

Carrier number of solitary laser A^soi 1.707 x 10^ 

Cavity decay rate F 0.357 x 10^^ s"^ 



n\t) = ]^G^n\t) [n\t)-ar{t)\ 

+ (Ji n\t - Ti) cos(a; Ti) + 
+ (Ji X'^(t — Ti) sin(cj Ti) + 

N 



CF2 ^ Gij W {t- T2) cos(cj T2) - 

N 

CF2 y^^GijV{t - T2) sin(cjr2) 



(74) 



V{t) = ]^G^n\t)[V{t)^an\t)\- 

— G\ VJ(t — Ti) sin(co' Ti) — 

— G\ X^(t — Ti) C0S(CJ Ti) — 

N 



- cr2 ^ Gij W {t- T2) sin(cj T2) - 

N 

- ^ Gij V{t - T2) cos(cj T2) (75) 

n\t) = {p-l)Jti,-jn\t)- 

-[T^G^n\t)]{[n\t)]'^[r{tf} (76) 

We integrate this differential equation system numeri- 
cally using Heun's method [56] which is a numerical inte- 
gration method of the class of Runge-Kutta methods that 
is particularly suitable for delay-differential equation sys- 
tems. The used step size for integration is At = 0.1 ps. In 
order to emulate the measurement of cross-correlations 
with a gigahertz oscilloscope in an experimental setup, 
on the one hand, we define cross-correlations between two 
simulated lasers as the cross-correlations between the ab- 
solute values of their electric fields. On the other hand, 
we use a sampling time of 1 ns per data point during 
which the absolute values of the electric field are aver- 
aged. For the precision measurement of the sensitivity 
of cross-correlations to detuning of the delay times in 
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Fig. 15, we use a sampling time of 10 ps. Desynchro- 
nization of the simulated lasers occurs during the power 
drop-outs of the low frequency fluctuations which hap- 
pen on a time scale of the order of magnitude of 10 r. In 
order to avoid measuring cross-correlations during them, 
the cross-correlations are calculated ten times from com- 
paratively short time windows of the length of 5r, and 
then the five largest cross-correlations are averaged. 

The master stability equations of the Lang-Kobayashi 
equations are 



1 



theorem defines determinants 



Sn{t) = -G^n{t) [Sn{t)-aSx{t)]- 



1 



Gn^uW [n{t)-al{t)]- 



+ CTi S-jzit - Ti) C0S(CJ Ti) + 

+ cTi 5x{t - Ti) sin(cjri) + 
+ cr2 Ik ^n{t - T2) cos(cj T2) + 
+ cr2 7fc ^x{t - T2) sin(a; T2) 

Sx{t) = J Gn n{t) [Sxit) + a Sn{t)] + 



(77) 



1 



cTi ^n{t - n) sin(a;ri) - 
(Ji Sx{t — Ti) cos(co'ri) — 
cr2 Ik ^n{t - T2) sin(a; T2) - 

Cr2 7/c ^x(t - T2) C0S(CJ T2) 



(78) 



2 [r + Gn n{t)] [n{t) 5n{t) + X{t) 5x{t)] (79) 



which are integrated numerically using Heun's method, 
as well. We calculate the maximal Lyapunov exponent 
using Farmer's method [40]. For correspondence to our 
definition of the cross-correlation above, we define the 
deviation of the absolute value of the electric field as the 
metric of the separation function. This is legitimate as 
the Lang-Kobayashi equations form a strongly coupled 
differential equation system. 



The Schur-Cohn theorem 

In order to find the synchronization areas in parameter 
space, the Schur-Cohn theorem is a possible method, see 
[38j.The synchronization area for a network of iterated 
one-dimensional maps with constant slope is derived by 
finding the parameters for which the roots of the charac- 
teristic polynomial lie inside the unit circle. 

For a polynom with P{x) = X]r=o Schur-Cohn 



an 








ao 


ai 


dn-l 


an 








ao 




dn-u+l ••• 


an 





.. 














oT 










■■ 











.. 



an 



(80) 



with i^ = 0,l,...,n — 1. 

Due to the Schur-Cohn theorem, it is |x| < 1 when all 
determinants are greater than 0. 
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